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ABSTRACT 

We present numerical solutions of the collapse of prestellar cores that lead to the formation and 
evolution of circumstellar disks. The disk evolution is then followed for up to three million years. 
A variety of models of different initial masses and rotation rates allows us to study disk accretion 
around brown dwarfs and low-mass T Tauri stars, with central object mass < 0.2 Mq, as well as 
intermediate and upper-mass T Tauri stars (0.2 < M* < 3.0 Mq). Our models include self-gravity 
and allow for nonaxisymmctric motions. In addition to the self-consistently generated gravitational 
torques, we introduce an effective turbulent a-viscosity with a = 0.01, which allows us particularly 
to model accretion in the low-mass regime where disk self-gravity is diminishing. A range of models 
with observationally-motivated values of the initial ratio of rotational to gravitational energy yields a 
correlation between mass accretion rate M and that is relatively steep, as observed. Additionally, 
our modeling reveals evidence for a bimodality in the M-M^ correlation, with a steeper slope at 
lower masses and a shallower slope at intermediate and upper masses, as also implied by observations. 
Furthermore, we show that the neglect of disk self-gravity leads to a much steeper M-M* relation for 
intermediate and upper-mass T Tauri stars. This demonstrates that an accurate treatment of global 
self-gravity is essential to understanding observations of circumstellar disks. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — ISM: clouds — stars: 
formation 



1. INTRODUCTION 

Using numerical hydrodynamic simulations, we 
have recently shown that self-gravitating disks of 
intermediate- and upper mass T Tauri stars (TTSs) set- 
tle into a self-regulated state, with low-amplitude den- 
sity perturbations persi sting for at least several Myr 
(|Vorobvov fc Basull2007t ). These perturbations are sus- 
tained by swing amplification at the disk's outer edge. 
The associated gravitational torques can drive mass 
accretion rates that arc of correct magnitude to ex- 
plain the observed values acr oss the range of interme - 
diate and upper mass TTSs (jVorobvov fc Basul [20081 ). 
These results were used to argue that the empirically 
observed correlation between mass accretion rate M 
and central object mass M*, of the approximate form 
M PC M? (e.g. iMuzerolle et al1l2003l. l2005l: iNatta et all 
[200l iCalvet et all 120041 : iMohantv et al.ll2005f ). can be 
explained purely on the basis of self-regulated accretion 
by gravitational torques in self-gravitating circumstellar 
disks. Nonaxisymmctric structure (and therefore nonax- 
isymmetric modeling) is required to achieve this result. 

However, it is important to note that the approxi- 
mate M oc empirical relation was only noticed after 
data were obtained for brown dwarfs (BDs) and com- 
bined with the TTS data. The larger dynamic range 
of masses and accretion rates in a combined plot re- 
veals a relatively steep overall gradient in M versus . 
Furthermore, there is evidence that the correlation is 
steeper at lower masses, and shallower at higher masses 
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(|Vorobvov fc Basull2008f ) , although this is complicated by 
the fact that different techniques are usually employed to 
determine M in the two m ass regimes (jMuzerolle et al.1 
[200l[Mohantv et al.l[200l . 

In this paper we model both the lower and upper mass 
end of the M-M^, correlation, using additional models 
for BD accretion that start from the collapse of very 
dense and compact low-mass cores. The low-mass disks 
that form around low-mass objects in our simulations are 
nearly axisymmetric and sustain negligible gravitational 
torques. Therefore, we explore here the possibility that 
additional processes are at play in the low-mass regime. 
These processes may depend on instabilities in the low- 
mass disks that generate turbulence that in turn yields 
an effective local viscosity. Our model is therefore con- 
sistent with the idea that gravitational torques are the 
driving agent of mass and angular momentum transport 
in intermediate- to high mass disks, while a mechanism 
like the magnetorotational instability (MRI) may dom- 
inate in low-mass disks. Other workers have arrived at 
this conclusi on based on physical and observational con - 
straints (e.g. iHartmann et al.ll2006l : iKratter et al]l2008[ ). 
We are also motivated by the apparent bimodal slope of 
the M-M* correlation, implying that two different phys- 
ical processes may be dominating disk evolution in the 
two different mass regimes. We note that our models are 
tracing the accretion processes for radii greater than 5 
AU, and we are implicitly assuming that the inner disk 
maintains the same accretion rate (at least in a time- 
averaged sense) through unspecified processes. There- 
fore, even the gravitational-torque-dominated disks may 
require processes such as the MRI to accomplish the cor- 
respondi ng angular momentum transport in the inner 
disk (see lArmitage et"alll2001l: [Zhu et al.ll2009[ ). 

We model the effect of turbulent viscosity with the a- 
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prescription (jShakura fc Sunvaevlll973l ). implemented in 
a nonaxisymmetric disk, as laid out in ^ It is applied to 
models of all masses. Protostellar disks are formed from 
the self-consistent collapse of initially slowly-rotating 
prestellar cores. The evolution of the disks are then fol- 
lowed for up to t hree million years. As explaine d in our 
previous papers (jVorobvov fc Basull2005l l2006f ). this is 
made possible through the use of the thin-disk approx- 
imation and our nonuniform grid. Such long-term cal- 
culations, including wide parameter surveys, remain out 
of reach for fully three-dimensional calculations. We use 
our results to compare the numerically-derived mass ac- 
cretion rates with those inferred observationally for TTSs 
and BDs in the late stage of star formation. 

2. MODEL DESCRIPTION 

2.1. Basic Equations 

Our numerical model is similar to that used recently 
to simulate the secular evolut ion of viscous and self - 
gravitating circumstellar disks (jVorobvov fc Bas"ijll2009f ) . 
Here, we briefly provide the basic concepts and equa- 
tions. We use the thin-disk approximation to compute 
the evolution of rotating, gravitationally bound cloud 
cores. The numerical integration is started in the prestel- 
lar phase, which is characterized by a collapsing starless 
cloud core, and continues into the late accretion phase, 
which is characterized by a protostar/disk/envelope sys- 
tem. This ensures a self- consistent formation of circum- 
stellar disks in our numerical simulations. Once the disk 
has formed, it occupies the innermost regions of our nu- 
merical grid, while the envelope occupies the rest of the 
grid. This means that the mass infall rate onto the disk 
is actually determined by the dynamics of gas in the en- 
velope, rather than being input through an ad hoc source 
term. The thin-disk approximation is an excellent means 
to calculate the evolution for many orbital periods and 
many model parameters. It is well justified as long as 
the ratio of the disk scale height Z to radius r does not 
considerably e xceed 0.1. As on e of us has recently shown 
(see figure 7 in I Vorobvovll2009l ) . this condition is fulfilled 
for solar mass stars having disks of several hundred AU 
in radius. Our model disks rarely exceed this size, and 
hence we believe that the thin-disk approximation is jus- 
tified after the disk formation epoch. This approximation 
is also reasonable for the prestellar p hase, since proto- 
stellar cores are found to be disk-hke (iJones et al.ll2001l : 
iJones fc Basull2002l: [Goodwin et al.ll2002l : lTassi£ll2007t ). 

The basic equations of mass and momentum transport 
in the thin-disk approximation are 

— = -V,.(S^,), (1) 
dv 

E^=-Vp7' + Sgp + (V.n)p, (2) 

where S is the mass surface density, V = jf^ Pdz is 
the vertically-integrated form of the gas pressure P, Z 
is the radially and azimuthally varying vertical scale 
height, Vp = VrV + v,f,(p is the velocity in the disk plane, 
= g^r + g^4> is the gravitational acceleration in the 

disk plane, and Vp = rd/dr + cj)r^^d/d<f> is the gradient 
along the planar coordinates of the disk. The gravita- 
tional acceleration includes both the gravity of a cen- 
tral point object (when formed) and the self-gravity of a 



circumstellar disk and envelope. The latter component 
is found by solving t he Poisson integral using the con- 
volution theorem (see IVorobvov fc Basul 120061 for more 
details). The viscous stress tensor is 

n=:2Si/^Vu-i(V-w)e^ , (3) 

where Vw is a symmetrized velocity gradient tensor, e is 
the unit tensor, and v is the kinematic viscosity. Equa- 
tion ([2]) describes the motion of a viscous fluid in the 
most general form. This equation can be reduced to the 
usual equation for the conservation of angular momen- 
tum of a radial annulus in the axisymmetric accretion 
disk. The component s of (V • II)p in p o lar co ordinates 
(r, (f)) can be found in IVorobvov fc Basul (|2009f ). 

It is well known that standard coUisional viscosity 
(molecular viscosity) is negligible in application to cir- 
cumstellar disks. An alternative is an effective turbu- 
lent viscosity induced by instabilities in the disk. A 
prime candidate is the MRI, although other mecha- 
nisms of nonlinear hydrodynamic turbulence cannot be 
completely ruled out due to the large Reynolds num- 
bers involved. We make no specific assumptions about 
the source of turbulence and parameterize the magni- 
tude of turbulent viscosity using the usual a-prescription 
(|Shakura fc SunvaevHTOTl 

V = aCsZ, (4) 

where = dV/dT, is the effective sound speed of (gen- 
erally) non- isothermal gas. The vertical scale height Z 
is determined in every computational cell and at every 
time step of integration using an assumption of local hy- 
drostatic equilibrium in th e gravitational field of th e cen- 
tral star and the disk (see IVorobvov fc Basull2009f ). It is 
thus important to note that our disks are not razor-thin 
but have a vertical extent and flaring which are consis- 
tent with those predicted from detailed vertical structure 
m odels of irradiated accre tion disks around T Tauri stars 
bv lD'Alessio et al.l l|1999( l. 

In this paper we present results for a large number of 
simulations with differing values of physical parameters 
but a fixed value a = 0.01. This choice is motivated 
by our recent analysis of t he secular evolution of v iscous 
and self-gravitating disks (|Vorobvov fc Basul |2009[ ). We 
found that i/ circumstellar disks around solar-mass proto- 
stars can generate and sustain turbulence then the tem- 
porally and spatially averaged a should lie in the range 
10"'^ — 10~^, in order to be a significant process but not 
completely destroy the self-regulated structure. There- 
fore, a choice of a = 0.01 samples models in which the 
turbulent viscosity has a significant and even dominant 
role during the late stages of disk accretion (gravitational 
torques are still dominant in the early stages for all but 
the lowest mass disks, as discussed in §|31). The possible 
influence of different values of a on our results is also 
discussed in Section [H 

Equations ([T|) and ^ are closed with a barotropic 
equation that makes a smooth transition from isother- 
mal to adiabatic evolution at S = !]„• = 36.2 g cm~^: 

V^clY. + clY.„[^ , (5) 

where Cs = 0.188 km s^^ is the sound speed in the 
initial state, and 7 = 1.4. Numerical simulations with 
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7 = 1.67 produce disks that are hotter than those with 

7 = 1.4. However, we have recently demonstrated that 
a modest increase in disk temperature has an insignifi- 
cant effect_on_the_jnBanjimss accretion rates (see figure 

8 in lVorobvov fc Basul [200i). Colder disks (T - 40 K at 
10 AU) are more gravitationally unstable and arc charac- 
terized by spiral modes that are of larger amplitude than 
those of hotter disks (T ~ 100 K at 10 AU). For hotter 
disks, however, the amplitudes of modes also decrease 
sharply with increasing order, leaving low-order modes 
(m < 2) dominant. As a result, the hotter models have 
less mode-to-mode interaction that can produce some 
cancellation in the net gravitational torque. Since the 
low-order modes are more efficient transport agents, the 
net effect is to produce comparable mass accretion rates 
in disks that differ in temperature by a factor of two. 
A similar result of increasing low-o rder mode domin ance 
for hotter disks has been found by ICai et all ()2008D . We 
further note that 7 depends on the disk temperature and 
hence it may vary radially, with the inner disk being char- 
acterized by 7 = 1.4 but the outer disk having 7 — 1.67. 
This radial variation in 7 is expected to decrease the ra- 
dial temperature gradient in the disk (as compared to the 
purely 7 = 1.4 disk) but is not expected to considerably 
change the mass accretion rates for the reason outlined 
above. 

2.2. Initial Conditions and Numerical Technique 

We start our numerical simulations from starless cloud 
cores, which have surface densities S and angular veloci- 
ties n typical for a collapsing axisymmetric magnetically 
supercritical core (|Basulll997D : 

ro^o 




(6) 



(7) 



where fio is the central angular velocity, ro is the radial 
scale length defined as ro = kcl/{GT,o) and k = \/2/tt. 
These initial profiles are characterized by the important 
dimcnsionless free parameter 7] = figrg/c^ and have the 
property that the asymptotic (r 3> ^o) ratio of centrifu- 
gal to grav itational acceleration has magnitude ^/2r] (see 
iBasulllQQTl ). The centrifugal radius of a mass shell that 
encloses a mass m and is initially located at radius r is 
estimated to be = p /{Gm) — \/2r\ r, where .7 = fir^ 
is the specific angular momentum (see lBasulll998D . We 
note that ry is similar in magnitude to the ratio of rota- 
tional to gravitational energy /? = i?rot/^'giav, where the 
rotational and gravitational energies are defined as 



-Ere 



27r 



ra 



-Itx 



rgr S r dr. 



(8) 

Here Uc = fl^r is the centrifugal acceleration, and rjn 
and rout are the inner and outer cloud core radii, respec- 
tively. The former is always set to 5 AU, while the latter 
varies according to the adopted cloud core size. From 
here onwards, we will refer to rj and /3 as synonymous 
quantities. The gas has a mean molecular mass 2.33 mn 



and cloud cores are initially isothermal with temperature 
T = 10 K. 

We present results from ten sets of models, each set 
being characterized by a distinct value of p. Individual 
models within every set are generated by varying rg and 
Hq in such a way that the product roilo is kept constant. 
This ensures that /3 is also constant for every model in the 
set, because the initial sound speed Cs = 0.188 km s~^ is 
equal for all models. All models have the ratio rout/''o 
set to 6.0 to generate truncated cloud cores of similar 
form. The values of /3, typical intervals for rout, ''0, ^0, 
and cloud core masses Md, and number of individual 
models within each set are listed in Table [1] We note 
th at our adopted value s of (3 lie within the limits inferred 
bv ICaselli et al.l (120021 ) for dense molecular cloud cores: 
/3 = (10-'' - 0.07). 

Equations dT]), ((21), ([5]) are solved in polar coordinates 
(r, 0) on a numerical grid with 128 x 128 points. We have 
found that an increase in the resolution to 256 x 256 grid 
zones has little influence on the accretion history. There- 
fore, we use the 128 x 128 grid in order to save a con- 
siderable amount of CPU time and to find solutions for 
a large number of model cloud cores. Each model takes 
about 400 CPU hours on the Opteron 2.5 GHz proces- 
sor. We use the method of finite differences with a time- 
explicit, operator-split solution procedure. Advection is 
performed using the second-order van Leer scheme. The 
radial points are logarithmically spaced. The innermost 
grid point is located at r^n ~ 5 AU, and the size of the 
first adjacent cell varies in the 0.17-0.36 AU range de- 
pending on cloud core size. We introduce a "sink cell" at 
r < 5 AU, which represents the central star plus some cir- 
cumstellar disk material, and impose a free inflow inner 
boundary condition. The outer boundary is reflecting. 
A small amount of artificial viscosity is added to the 
code, though the associated artificial viscosity torques 
were shown to be negligible in compari son with gravita- 
tional torques (jVorobvov fc Basull2007[ ). For reasonable 
values of the a-parameter (> lO"'' — 10"'^), the artificial 
viscosity torques are also considerably smaller than the 
viscous torques due to a- viscosity. 

3. RESULTS 
3.1. Time Evolution 

We begin with reviewing the accretion history in ob- 
jects formed from cloud cores of distinct masses but 
having equal ratios of the rotational-to-gravitational en- 
ergy p. For this purpose, we choose model set 6 with 
P = 8.0 X 10"'^. The corresponding mass accretion rates 
versus time elapsed since the beginning of simulations 
are shown in two upper rows of Figure [T] The time- 
dependent mass accretion rate is calculated as M{t) = 
— 27rrinWrS, where Vr is the inflow velocity through the 
sink cell and r-^ = 5 AU is the radius of the sink cell. The 
initial cloud core mass Md and central angular velocity 
fig are shown in each frame. 

The early accretion history for all objects is similar — 
M reaches a peak value of ~ 2.0 x 10^^ Mq yr^^ soon 
after the central object formation and settles at a near 
constant value of ~ 10^^ Mq yr^^. Then, a transient 
sharp decline in M follows, manifesting the onset of disk 
formation. For a short time, centrifugal forces balance 
those of viscosity and gravity and the mass accretion 
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TABLE 1 

Model parameters 



Set 


/3 




no 








rout 


Mci 


N 


1 


6.2 X 10~ 




0.28 


- 0.47 


2070 


- 3450 


(1.2 


- 2.0) X 10* 


1.2 - 


2.0 


5 


2 


1.3 X IQ- 


3 


0.27 


- 1.33 


1040 


- 5180 


(0.6 


- 3.0) X lO"' 


0.6 - 


3.0 


11 


3 


2.5 X IQ- 


3 


0.37 


- 3.40 


570 


- 5180 


(0.3 


- 3.0) X 10-1 


0.3 - 


3.0 


16 


4 


3.7 X IQ- 


3 


0.52 


- 6.80 


350 


- 4490 


(0.2 


- 2.6) X lO** 


0.2 - 


2.6 


18 


5 


5.1 X IQ- 


3 


0.53 


- 10.0 


280 


- 5180 


(0.16 


- 3.0) X 10"' 


0.16 - 


- 3.0 


17 


6 


8.0 X IQ- 


3 


0.67 


- 20.0 


172 


- 5180 


(0.1 


- 3.0) X 10" 


0.1 - 


3.0 


18 


7 


1.2 X IQ- 


2 


0.80 


- 34.0 


120 


- 5180 


(0.07 


- 3.0) X lO" 


0.07- 


- 3.0 


16 


8 


1.6 X IQ- 


2 


2.30 


- 50.0 


95 - 


- 2070 


(0.06 


- 1.2) X 10" 


0.05 - 


- 1.2 


12 


9 


2.0 X IQ- 


2 


1.60 


- 91.0 


60 - 


- 3450 


(0.04 


- 2.0) X 10" 


0.03 - 


- 2.0 


14 


10 


3.2 X IQ- 


2 


2.0 - 


- 114.0 


60 - 


- 3450 


(0.04 


- 2.0) X 10" 


0.03 - 


- 2.0 


13 



Note. — All distances are in AU, angular velocities in km s ^ pc 

rate drops to a negligible value. As the disk contin- 
ues to build up its mass, the burst mode of accretion 
ensues. During this phase, prolonged periods of rela- 
tively low accretion at 10~^-10~^ Mq yr~^ are inter- 
spersed with short episodes of activity when M increases 
to 10~^-10~^ Mq yr^^. These accretion bursts are as- 
sociated with disk fragmentation and formation of dense 
massive clumps that are later driven onto the central 
star due to the gravitatio nal interaction with spiral arms 
(jVorobvov fc Basul I2006D . Figure [T] demonstrates that 
the intensity and duration of the burst mode increases 
along the line of increasing cloud core masses. Most of 
the burst activity is constrained to the early several hun- 
dred thousand years of evolution, though the most mas- 
sive objects can undergo a few bursts even after 0.5 Myr. 
In this late evolution phase, the mass accretion rates 
show a moderate decline by roughly one order of magni- 
tude during 2-3 Myr. Some short-term variability is also 
present, but its magnitude is much smaller than that of 
the long-term decline. 

Mass accretion rates in other model sets show a simi- 
lar pattern of behavior. We now turn to comparing our 
numerically derived accretion rates in the late evolution 
phase with those inferred from observations. 

3.2. Mass Accretion Rates 

In this section, we compare our model mass accretion 
rates and central obje ct masses with those compiled by 
iMuzeroUe et all (|2005l . and references therein) for TTSs 
and BDs of age 0.5-3.0 Myr based on the measurements 
in (mostly) Taurus and other star formation regions. In 
order to facilitate the comparison, we time-average our 
model mass accretion rates and central object masses 
between 0.5 Myr and 3.0 Myr from the central object 
formation to obtain characteristic mean values (M) and 
(M*), respectively. For each model in Figure [Tl we 
also construct the distribution function (DF) of accre- 
tion rates by calculating M every 20 yr between 0.5 Myr 
and 3.0 Myr from the central object formation and dis- 
tributing the resulted values among 100 logarithmically- 
spaced bins in the 10~^ M© yr~^-10~^ Mq yr~^ range. 
The resultant DFs representing the normalized number 
of accretion "measurements" N-^ (by analogy to obser- 
vations) in a given mass accretion bin are shown in two 
bottom rows of Figure [TJ It is seen that the accretion 
rates in each model are localized in a specific interval. 
There is a preference for lower values of M in each inter- 
val, implying the measurements of mass accretion rates 



^ , and masses in Mq . N is the number of models in each model set. 

may be biased towards lower values. 

Color symbols in Figure [5] represent the time-averaged 
mass accretion rates (M) (in Mq yr~^) versus time- 
averaged object masses (M*) (in Mq). Each filled symbol 
(of same color and shape) within a given set of models 
represents an individual object, which has formed from 
a cloud core of distinct mass, angular velocity, and size. 
In addition, each filled symbol is assigned vertical bars 
showing typical variations in the accretion rate between 
0.5 Myr and 3.0 Myr after the central object formation. 
The typical variations are defined as mass accretion rates 
with Ny^ equal to or greater than 0.05. This procedure 
helps to exclude objects with accretion rates that have a 
probability to be detected that is less than one part in 
twenty and are thus statistically insignificant. Consid- 
ering the number (< 20) of detected objects along the 
line of constant central object mass, we believe that our 
typical variations represent conservative estimates. The 
open symbols and plus signs in Figure [2] represent data 
obtained from observations. In particular, the open di- 
amonds represent meas urements, mostly in Tau rus, that 
have been compiled by iMuzeroUe et al.l (j2005l . and ref- 
erences therein); the open square s repre sent detections 
in p Oph obtained bv lNatta et al.l ()2006f ) and plus signs 
represent their upper limits to nondete ctions. We have 
exclud ed objects in the compilati on of | Muzerolle et al.l 
()2005f) that were later observed bv lNatta et al.l 1)20061 ). 

There are several important conclusions that can be 
drawn by analyzing Figure [H 

1. Our numerical models yield mass accretion rates 
that are of correct magnitude to explain the ob- 
served values in TTSs and BDs. There are only a 
few objects whose accretion rates fall beyond our 
predicted limits, but these objects can easily be ac- 
counted for by increasing /3 marginally above our 
adopted values. 

2. The time- averaged mass accretion rates and stellar 
masses for objects from every individual set of mod- 
els form a unique track. These tracks are distinct 
for BDs and low-mass TTSs but tend to converge 
for upper-mass TTSs. Objects within each individ- 
ual track have formed from cloud cores of distinct 
mass. In particular, objects in the lower-left end 
of the track have formed from cloud cores of lower 
mass. 

3. Our models suggest that TTSs of 0.5-3.0 Myr age 
may have a wider range of mass accretion rates 
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Fig. 1. — Two upper panels. Mass accretion rates versus time elapsed since the beginning of simulations. Six models from model 
set 6 arc presented, the initial cloud core mass M^i (in Mq) and central angular velocity Co (in k™ s~^ pc~^) are indicated in every 
panel. Vertical dotted lines mark the onset of the late accretion phase at 0.5 Myr from the central object formation. Two lower 
panels. Theoretical distribution functions (DFs) of accretion rates for the same six models showing the normalized number of accretion 
"measurements" as a function of the mass accretion rate. The measurements are done every 20 yr between 0.5 Myr and 3.0 Myr from the 
central object formation. 



than implied from observations. In particular, we 
predict the existence of TTSs with accretion rates a 
factor of 10 lower than has been fo und in the com- 
pilation of iMuzeroUe et al.1 (|2005l . and references 
therein). The lack of such objects in the observa- 
tional data is most likely explained by a difficulty 
to measure extremely low accretion rates in TTSs. 

A dearth of BDs with accretion rates M > 
10~^ Mq yr~^ is likely caused by the lack of cloud 
cores with sufficiently large rotation rates. For in- 
stance, in order to reproduce the largest detected 
M in BDs (« 10~* Mq yr'^), we had to employ 
models with /3 = 0.032. This ratio of the rotational 



to gra vitational energy, according to ICaselli et al.l 
()2002f ). is close to the upper measured limit for 
dense cloud cores. 

Perhaps, the most important conclusion from our mod- 
eling is that the observed scatter in the mass accretion 
rates along the line of equal object masses cannot be ex- 
plained by intrinsic variability (shown in Figure[2]by ver- 
tical bars). This variability can account only for a max- 
imum of two orders of magnitude scatter (with the no- 
ticeable exception of a few low-mass BDs) and often even 
less, but the measured rates typically span a range of 
three orders of magnitude. Figure [5] clearly demonstrates 
that some object-to-object variations are necessary to ex- 
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0.1 1 

Central object mass (Mq) 

Fig. 2. — M a ss acc retion rates versus central object masses. The open diamonds represent measurements of TTSs and BDs from 
IMuzeroUe et al] II2005I . and references therein); the open squares and plus signs represent the confirmed detections and upper limits, 
respectively, compiled bv INatta et al.l II2006I ). The filled color symbols show the time-averaged mass accretion rates (M) versus time- 
averaged central object masses (M*) obtained in our numerical modeling. Each symbol type represents a distinct set of models with a 
constant value of /3, as detailed in Table [T] In particular, model set 1 is plotted by cyan diamonds, model set 2 — by red trianglcs-up, 
model set 3 — by black circles, model set 4 — by green squares, model set 5 — by red circles, model set 6 — by blue triangles-down, model set 
7 — by pink circles, model set 8 — by cyan circles, model set 9 — by yellow circles, and model set 10 — by red diamonds. The bars represent 
typical variations in accretion rates in each model. The solid/dashed lines are the least-squares best fits to the model/observational data, 
respectively. In particular, the left lines are the fits to BDs and low-mass TTSs and right lines are the fits to the intermediate- and 
upper-mass TTSs. 



pl ain the observed sca tter. This confirms previo us claims 
bv INatta et all (|2004f ) and lNguven et all ()2009f ). 

The object-to-object variations in the mass accretion 
rate along the line of equal (sub)stellar masses are caused 
by the difference in the disk masses. More massive 
disks are expected to d rive higher rates of accretion 
() Vorobvov fc Basul [20081 ) . Figure H shows that objects 
with greater (M) but equal (A/*) belong to model sets 
characterized by greater values of /3. At the same time, 
P controls the disk mass because the centrifugal radius 
is directly proportional to (see Section 12. 2p . This 
means that models with greater /3 (and greater (M)) 
but equal (M*) are generally expected to harbor more 
massive disks and this is c onfirmed by our calculations 
of disk masses presented in IVorobvovl ()2009[ ) . 

When all observational data in Figure [2] for TTSs and 
BDs are taken together, a least squares fit is described 
by a power law 

M = 10-^-X'-'^°-'- (9) 

When wc take the least-squares fit to the confirmed de- 
tections only (excluding non-detections represented by 
plus signs), we obtain the relation 

M = 10-^-5m2°±oi. (10) 

The difference between the best fits is not significant, 
but it demonstrates that the exponent is sensitive to the 



way the data are handled and may actually change when 
more observational data become available. 

The least-squares best fit to all model data shown in 
Figure [2] is described by the relation 

(M) = 10-^-3(M,)i-8±°-i. (11) 
It is evident that our numerical modeling can reproduce 
the observed relation reasonably well, except probably 
for the fact that we somewhat overestimate the observed 
rates. Possible reasons for this small discrepancy are 
discussed in Section 

However, we believe that taking a best fit over the 
whole mass range of BDs an d TTSs may be mislead - 
ing. In our previous paper l| Vorobvov fc Basul |2008[) . 
we reanalyzed the data obtained from observations and 
argued that the exponent n in the M oc relation 
takes different values for the lower-mass and upper-mass 
objects. In particular, we found n = 2.3 ± 0.6 for 
BDs and low-mass TTSs with < 0.25 Mq, and 
n = 1.3 ± 0.3 for intermediate- and upper mass TTSs 
with 0.25 Mq < < 3.0 Mq. Two possible explana- 
tions were put forward for this apparent bimodality: vari- 
ations in the observational methods and different mech- 
anisms responsible for accretion. In this paper we focus 
on the latter alternative. 

The left and right solid lines in Figure H] arc the least- 
squares fits to our model data for the lower- and upper- 
mass objects, respectively. The dashed lines are the cor- 
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responding best fits to the observational data. The least- 
squares fits in each mass regime are distinct and are de- 
scribed by the following power laws. 
Modeling: 

r (A/) = fO-6-3 (Af,)2-9±0-5 if (M,) < 0.2 Mq, 
{ (A/) = f 0-^-3 (A4)i-5±0-i if 0.2 Mq < < 3.0 

(12) 

Observations: 

r M = 10-7-2 7Vf2.3±0.6 if < 0.2 Mq, 

\ M ^ 10-^-6 Afi-3±0-3 if 0.2 Mq < M^ < 3.0 Mq . 

It is evident that our numerical modeling corroborates 
the presence of bimodality in the observed M-M^, re- 
lation, though predicting a steeper dependence of M on 
Af* for BDs and low-mass TTSs. We note, however, that 
a shallower dependence in the low-mass regime could be 
obtained by expanding our modeled range of values of /3. 

4. DISCUSSION 

How can different values of a affect our model accretion 
rates shown in Figures [T] and [2]? Let us first focus on 
the intermediate- and upper mass TTSs with 0.2 Mq < 
M^ < 3.0 Mq. T Tauri disks are unlikely to sustain a > 
0.1 throughout the whole disk volume during a typical disk 
lifetime of 3-5 Myr. Such large values of a lead to a disk 
deple tion and disappearance on time scales of less than 1 
Myr (jVorobvov fc Basull2009l ) . The results of that study 
excludes large values of a and leaves us with a < 10"^ as 
the most probable range of values. T he accretion rates 
in the a = case were presented by iVorobvov fc Basul 
(|2008[) . The (M) oc (Af*)" relation in this case was found 
to have an exponent n ~ 1.7 ± 0.1, which is similar to 
the value found in the present study (see Equation [T^. 
but the values of (A/) were approximately a factor of 2 
smaller than in the a = 10^^ case. This means that the 
actual value of a will not dramatically change the mean 
mass accretion rates, as long as it is < 10"^. 

The effect of varying a in the BD and low-mass TTS 
regime is more difficult to assess. Circumstellar disks in 
this mass regime have a lower surface density (than disks 
around intermediate- and upper-mass TTSs) and hence 
are expected to be more MRTactive. On the other hand, 
we do not expect a to be equal to or greater than 0.1 due 
to the same reasons as discussed above. The fact that 
the least-squares fits to the model and observational data 
in Figure [2] tend to converge at the low-mass end implies 
that a = 0.01 is likely to be a good choice, certainly for 
the BD mass regime. The limit of small a < 10""^ is dif- 
ficult to examine due to numerical reasons. Normally, as 
in the case of intermediate- and upper-mass TTSs, disks 
are relatively massive and artificial vi scosity torques are 
much smaller than gravitational ones (jVorobvov fc Basul 
|2007| ). However, when both a and disk masses are small, 
artificial viscosity torques may become comparable to 
both gravitational and viscous ones. In addition, the 
imperfections at the inner inflow boundary could also in- 
troduce some low amplitude noise in the inner regions, 
which may not be negligible in the limit of small a and 
small disk masses. In this situation, the model mass ac- 
cretion rates may be artificially overestimated. To avoid 
these complications, we decided to use a = 10~^ in our 
modeling. 



Central object mass (M,^,) 
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open squares - viscosity plus disk self-gravity ! 
filled squares - viscosity only 
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Central object mass (Mg.) 

Fig. 3. — Time-avcragcd mass accretion rates (Af ) versus time- 
averaged central object masses (Af, ) for all models from model set 
3 (top) and model set 6 (bottom). Open squares show the case 
when both disk self-gravity and viscosity are at work, while filled 
squares correspond to the case when disk self-gravity is set to zero 
and viscosity remains the only mass transport agent in the disk. 
The solid and dashed lines are the least-squares best fits to models 
with disk self-gravity and without disk self-gravity, respectively. 

In the early disk evolution at t < 0.5 Myr from the 
central star formation, gravitational-instability-induced 
torques dominate viscous ones, for interrnediate and 
upper-mass TTSs (see e.g. IVorobvov fc Basu|[2009t ). and 
this can influence considerably the subsequent time be- 
havior of the system. We illustrate this phenomenon in 
Figure|31 which shows time-averaged mass accretion rates 
(M) versus time-averaged central object masses (Af*) for 
model set 3 (top) and model set 6 (bottom). Every model 
in these sets is run twice: first time — with both disk self- 
gravity and viscosity and second time — with disk vis- 
cosity only. In the latter case, the disk self-gravity is 
artificially set to zero, thus nullifying the gravitational 
instability as well, but the gravity of the central object 
is kept. The open squares in Figure |3| show the data ob- 
tained in models with both disk self-gravity and viscos- 
ity (these data are identical to those shown in Figure 
while the filled squares present the data for models with 
disk viscosity only. 

It is evident that gravitational instability has little ef- 
fect on the low-mass objects in each model set, but its 
effect becomes noticeable along the line of increasing 
object masses. Starting from the fourth or fifth least 
massive object in each set, the models without disk self- 
gravity begin to yield an (M) that is greater than those 
of self-gravitating models, while at the same time un- 
derestimating the central object mass (Af*). In general, 
models without self-gravity tend to produce stars with 
masses that arc too low and mass accretion rates that 
are too high. This effect may seem counterintuitive since 
both viscous and gravitational torques are expected to 
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act together to increase {M). However, it is important 
to keep in mind that our time-averaged mass accretion 
rates apply to the late evolution phase between 0.5 Myr 
and 3.0 Myr after the formation of the protostar. In 
contrast, during the early phase (< 0.5 Myr), the time- 
averaged mass accretion rates in models with self-gravity 
are greater than in models without self-gravity. This is 
due to vigorous gravitational instability that acts in self- 
gravitating disks. As a result, models with self-gravity 
build up the stellar mass and deplete the total gas reser- 
voir much faster than models without self-gravity. This 
causes the mass accretion rate in self-gravitating disks to 
be lesser in the late evolution phase. 

In order to better illustrate the difference in the mass 
accretion rates, we take the least-squares best fits to the 
models with self-gravity (solid lines) and models with- 
out self-gravity (dashed lines). In the former case, we 
obtain exponents n = 2.5 ± 0.2 (top panel, model set 3) 
and n = 2.1 ± 0.2 (bottom panel, model set 6). When 
the least-squares fit is taken for both sets of models, we 
obtain n = 2.0 ± 0.2. This is quite similar to the expo- 
nent found for the complete set of models (see eq. [H]). 
In the case without self-gravity, we obtain exponents 
n — 4.2 ± 0.2 (top panel, model set 3) and n = 3.6 ± 0.2 
(bottom panel, model set 6), which are considerably 
larger than those for models with self-gravity. In this 
example, we have not considered all models due to an 
enormous computational load, but we do not expect this 
tendency to change considerably for the complete set of 
models. This allows us to conclude that if not for the disk 
self-gravity and associated gravitational instability in the 
early evolution phase, we would have had difficulty to 
recover the observed M~M^ relation. 

Another important feature of the non-self-gravitating 
models in Figure [3] is an apparent lack of bimodality in 
the M-M^ relation. Almost all models within an individ- 
ual set, except for a few ones with the smallest masses 
of the central object, fall onto a track that is well de- 
scribed by a straight line in log-log space. This strongly 
suggests that the bimodality is the result of gravitational 
instability acting in the early phase of disk evolution. 

5. CONCLUSIONS 

Using numerical hydrodynamic simulations of circum- 
stellar disk formation and evolution, we have studied 
the mass accretion rates in BDs and TTSs, focusing 
mainly on the late evolution phase between 0.5 Myr and 
3.0 Myr from the central object formation. Our numer- 
ical model involves both disk self-gravity, which is accu- 
rately computed via the solution of the Poisson integral, 
and turbulent viscosity, the latter being described by a 
usual Shakura & Sunyaev parameterization with the a- 
parameter set to 10"'^. The theoretical data are com- 
pared against those obtaine d from the measuremen ts in 
young star-forming clusters iMuzeroUe et al.l ()2005l . and 
references therein). We find the following. 

• Our numerical modeling of self-consistently-formed 
circumstellar disks yields mass accretion rates of 
correct magnitude to explain the observed values 
in TTSs and BDs of 0.5-3.0 Myr age. 



• We corroborate our previous conclusion 
(jVorobvov fc Basul l2008f ) that the dependence 
of the mass accretion rates (M) on the central 
object mass (M*) in TTSs and BDs can be 
better described by a bimodal power-law func- 
tion rather than that with a single exponent. 
Mass accretion rates of BDs and low-mass TTSs 
(M* < 0.2 Mq) have a steeper dependence on Af* 
than those of the intermediate- and upper-mass 
TTSs (0.2 Mq < < 3.0 Mq). In particular, the 
least-squares fits to our model data yield exponents 
n = 2.9 ±0.5 and n — 1.5 ±0.1 for the objects with 
mass Mh. < 0.2 Mq and 0.2 Mq < < 3.0 Mq, 
respectively. The corresponding fits to the obser- 
vational data produce exponents n ~ 2.3 ± 0.6 and 
n = 1.3 ±0.3. 

• The apparent bimodality in the M-M^ relation is 
caused by vigorous gravitational instability in the 
early phase of disk evolution. The gravitational in- 
stability serves to limit disk i nasses in the int erme- 
diate and upper-mass TTSs (|Vorobvovll2009l) . thus 
effectively setting an upper limit on the mass accre- 
tion rates in the late evolution and flattening the 
M-M^ relation in this mass regime. Models with- 
out self-gravity greatly overestimate the observed 
mass accretion rates, while at the same time under- 
estimating the central object masses. As a result, 
the non-self-gravitating models fail to account for 
the observed M-M* relation, predicting a much 
steeper dependence of M on M*. 

• The observed large scatter in mass accretion rates 
along the line of equal (sub)stellar masses is caused 
in part by the intrinsic variability during the evolu- 
tion of individual objects, which may span a range 
of one (intermediate- and upper-mass TTSs) to two 
(BDs and low-mass TTSs) orders of magnitude. 
The other part is object-to-object variations due 
to different initial conditi ons. Our conclusi on is 
in agreement with t hat of iNatta et al.l (|2004f ) and 
iNguven et al.l (j2009f ) based on the analysis of vari- 
ability in TTSs. 

• We predict the existence of TTSs with accretion 
rat es a factor of 1 low er than those reported 
bv iMuzeroUe et all (j2005l . and references therein). 
Our modeling also indicates that upper mass TTSs 
may exhibit a few FU-Ori-like outbursts with ac- 
cretion rates greater than 10~^ Mq yr~^. However, 
the anticipated number of such peculiar objects is 
quite small and they are not expected to break the 
existing M-M* relation. 



The authors are thankful to the anonymous referee for 
very helpful suggestions. EIV gratefully acknowledges 
support from an ACEnet Fellowship. SB was supported 
by a grant from NSERC. We thank the Atlantic Compu- 
tational Excellence Network (ACEnet) and the SHARC- 
NET consortium for access to computational facilities. 



REFERENCES 

Armitage, P. J., Livio, M., & Pringlc, J. E. 2001, MNRAS, 324, Basu, S. 1997, ApJ, 485, 240 

705 



T Tauri and Brown Dwarf Disk Accretion 



9 



Basu, S. 1998, ApJ, 509, 229 

Cai, K., Duriscn, R. H., Boley, A. C, Pickett, M. K., & Meji'a, A. 

C. 2008, ApJ, 673, 1138 
Calvet, N., MuzeroUe, J., Briceno, C, Hernandez, J., Hartmann, 

L., Sauceno, J. L., & Gordon, K. D. 2004, AJ, 128, 1294 
Caselli, P., Benson, P. J., Myers, P. C, & Tafalla, M. 2002, ApJ, 

572, 238 

D'Alessio, P., Calvet, N., Hartmann, L., Lizano, S., & Canto, J. 

1999, ApJ, 527, 893 
Goodwin, S. P., Ward-Thompson, D., & Whitworth, A. P. 2002, 

MNRAS, 330, 769 
Hartmann, L., D'Alessio, P., Calvet, N., & MuzeroUe, J. 2006, ApJ, 

648, 484 

Jones, C. E., Basu, S., & Dubinski, J. 2001, ApJ, 551, 387 
Jones, C. E., & Basu. S. 2002, ApJ, 569, 280 

Kratter, K. M., Matzner, C. D., & Krumholz, M. R. 2008, ApJ, 
681, 375 

Mohanty, S., Jayawardhana, R., & Basri, G. 2005, ApJ, 626, 498 
MuzeroUe, J., Hillenbrand, L., Calvet, N., Bricerio, C, & 
Hartmann, L. 2003, ApJ, 592, 266 



MuzeroUe, J., Luhman, K., Bricerio, C, Hartmann, L., & Calvet, 

N. 2005, ApJ, 625, 906 
Natta, A., Testi, L., MuzeroUe, J., Randich, S., Comeron, F., & 

Persi, P. 2004, A&A, 424, 603 
Natta, A., Testi, L., & Randich, S. 2006, A&A, 452, 245 
Nguyen, D. C, Scholz, A., van Kerkwijk, M. H., Jayawardhana, 

R., Brandeker, A. 2009, ApJ, 694, L153 
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Tassis, K. 2007, MNRAS, 379, L50 
Vorobyov, E. I. 2009, ApJ, 692, 1609 
Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137 
Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956 
Vorobyov, E. I., & Basu, S. 2007, MNRAS, 381, 1009 
Vorobyov, E. I., & Basu, S. 2008, ApJ, 676, 139 
Vorobyov, E. I., & Basu, S. 2009, MNRAS, 393, 822 
Zhu, Z., Hartmann, L. & Gammic, C. F. 2009, ApJ, 694, 1045 



